Introducción

Para este trabajo de clasificación binaria utilizando diversas técnicas de machine learning, se ha decidido utilizar el set de datos Dengue prevalence by administrative region.

Este dataset contiene información relevante sobre la prevalencia de la enfermedad por cada región administrativa de Nueva Zelanda (2000 regiones), es decir, si se ha observado la presencia del dengue o no, desde 1961 hasta 1990.

El objetivo de este trabajo es entrenar diversos algoritmos de clasificación, con la finalidad de predecir una variable binaria y de tomar una decisión sobre cual modelo propuesto es el mas recomendado para el tipo de dataset con el que se está trabajando.

Fuentes

  1. Repositorio: https://vincentarelbundock.github.io/Rdatasets/datasets.html

  2. Explicación de los datos (en ingles): https://vincentarelbundock.github.io/Rdatasets/doc/DAAG/dengue.html

  3. Descripción de las columnas que contiene el dataset dengue:

    3.1 humid: Densidad de vapor promedio desde 1961 a 1990.

    3.2 humid90: Percentil 90 para humedad.

    3.3 temp: Temperatura promedio desde 1961 a 1990.

    3.4 temp90: Percentil 90 para temperatura.

    3.5 h10pix: Humedad máxima dentro de un radio de 10 pixeles.

    3.6 h10pix90: Humedad máxima de variable temp90 dentro de un radio de 10 pixeles.

    3.7 trees: Porcentaje de un área cubierta por arboles, dato entregado por satélite.

    3.8 trees90: Percentil 90 de la variable trees.

    3.9 NoYes: Se ha observado la presencia de dengue, 1 indica si. Variable dependiente a estudiar.

    3.10 Xmin: Longitud mínima.

    3.11 Xmax: Longitud máxima.

    3.12 Ymin: Latitud mínima.

    3.13 Ymax: Latitud máxima.

    3.14 X: Indicador fila.

  4. Dimension de los datos: 2000 filas x 14 columnas

Analisis exploratorio de los datos

X humid humid90 temp temp90 h10pix h10pix90 trees trees90 NoYes Xmin Xmax Ymin Ymax
1 0.6713889 4.416667 2.037500 8.470835 17.35653 17.80861 0 1.5 0 70.5 74.5 38.0 35.5
2 7.6483340 8.167499 12.325000 14.925000 10.98361 11.69167 0 1.0 0 62.5 64.5 35.5 34.5
3 6.9790556 9.563057 6.925000 14.591660 17.50833 17.62528 0 1.2 0 68.5 69.5 36.0 35.0
4 1.1104163 1.825361 4.641665 6.046669 17.41764 17.51694 0 0.6 0 67.0 68.0 35.0 34.0
5 9.0270555 9.742751 18.175000 19.710000 13.84306 13.84306 0 0.0 0 61.0 64.5 33.5 32.0
6 8.9141113 9.516778 11.900000 16.643341 11.69167 11.69167 0 0.2 0 64.5 65.5 36.5 35.0
'data.frame':   2000 obs. of  14 variables:
 $ X       : int  1 2 3 4 5 6 7 8 9 10 ...
 $ humid   : num  0.671 7.648 6.979 1.11 9.027 ...
 $ humid90 : num  4.42 8.17 9.56 1.83 9.74 ...
 $ temp    : num  2.04 12.32 6.93 4.64 18.18 ...
 $ temp90  : num  8.47 14.93 14.59 6.05 19.71 ...
 $ h10pix  : num  17.4 11 17.5 17.4 13.8 ...
 $ h10pix90: num  17.8 11.7 17.6 17.5 13.8 ...
 $ trees   : num  0 0 0 0 0 0 0 0 0 0 ...
 $ trees90 : num  1.5 1 1.2 0.6 0 ...
 $ NoYes   : int  0 0 0 0 0 0 0 0 0 0 ...
 $ Xmin    : num  70.5 62.5 68.5 67 61 64.5 67.5 64 63.5 61 ...
 $ Xmax    : num  74.5 64.5 69.5 68 64.5 65.5 68.5 66.5 65.5 64.5 ...
 $ Ymin    : num  38 35.5 36 35 33.5 36.5 33.5 35 33 35 ...
 $ Ymax    : num  35.5 34.5 35 34 32 35 32 33 29.5 33.5 ...
[1] 2000   14

Como se observa anteriormente, son 14 columnas con 2000 filas, los tipos de datos son de tipo numérico incluyendo la variable objetivo NoYes, se realiza a continuacion una exploración mas detallada del dataset a estudiar.

En cuanto a las observaciones sobre la variable independiente (NoYes) y observando el grafico de referencia, se observa un dataset un tanto desbalanceado pero no deberia afectar en demasia la ejecucion de esta practica.

NoYes cuenta porcentaje
0 1169 58.45
1 831 41.55

A continuación se presentan algunos gráficos de apoyo al análisis exploratorio:

Observaciones:

Transformaciones

Se realizan transformaciones sobre los datos, se eliminan columnas, tratamiento de datos faltantes y discretización de los datos:

# Se elimina la columna X
datos$X <- NULL

# Tratamiento de datos faltantes, se eliminan
datos2 <- na.omit(datos, (!is.na(datos)))
colSums(is.na(datos2)) > 0 
   humid  humid90     temp   temp90   h10pix h10pix90    trees  trees90 
   FALSE    FALSE    FALSE    FALSE    FALSE    FALSE    FALSE    FALSE 
   NoYes     Xmin     Xmax     Ymin     Ymax 
   FALSE    FALSE    FALSE    FALSE    FALSE 
# Se separa la variable objetivo de la variable input
varObjBin <-datos2$NoYes
input <- as.data.frame(datos2[, -(9)])
# Pre procesado, recategorizar, dummies, estandarizar etc

# Estandarizacion
temp <- input %>% mutate_if(is.numeric, scale)
head(temp)
      humid   humid90        temp      temp90     h10pix   h10pix90      trees
1 -2.187817 -1.752592 -2.04603997 -1.42684042 -0.5252414 -0.5182986 -0.9516725
2 -1.235570 -1.240199 -0.76267048 -0.58677814 -1.3943124 -1.3621731 -0.9516725
3 -1.326917 -1.049555 -1.43632247 -0.63016510 -0.5045405 -0.5435904 -0.9516725
4 -2.127896 -2.106585 -1.72116940 -1.74236542 -0.5169085 -0.5585368 -0.9516725
5 -1.047396 -1.025007 -0.03288072  0.03602869 -1.0043711 -1.0653740 -0.9516725
6 -1.062811 -1.055877 -0.81568941 -0.36312211 -1.2977547 -1.3621731 -0.9516725
    trees90      Xmin      Xmax      Ymin      Ymax
1 -1.173965 0.9177821 0.9477775 0.7759536 0.7328488
2 -1.191357 0.7886923 0.7859356 0.6691044 0.6904187
3 -1.184400 0.8855096 0.8668566 0.6904742 0.7116337
4 -1.205271 0.8613053 0.8425803 0.6477345 0.6692036
5 -1.226141 0.7644880 0.7859356 0.5836250 0.5843434
6 -1.219184 0.8209648 0.8021198 0.7118441 0.7116337
# Se une la variable objetivo con los resultados(de estandarizar, dummys etc)
base <- data.frame(varObjBin, temp)

# La variable objetivo se cambia a alfanumerico yes, no
base$varObjBin <- ifelse(base$varObjBin == 1, "Yes", "No")
knitr::kable(head(base), "pipe")
varObjBin humid humid90 temp temp90 h10pix h10pix90 trees trees90 Xmin Xmax Ymin Ymax
No -2.187817 -1.752592 -2.04603997 -1.42684042 -0.5252414 -0.5182986 -0.9516725 -1.173965 0.9177821 0.9477775 0.7759536 0.7328488
No -1.235570 -1.240199 -0.76267048 -0.58677814 -1.3943124 -1.3621731 -0.9516725 -1.191357 0.7886923 0.7859356 0.6691044 0.6904187
No -1.326917 -1.049555 -1.43632247 -0.63016510 -0.5045405 -0.5435904 -0.9516725 -1.184400 0.8855096 0.8668566 0.6904742 0.7116337
No -2.127896 -2.106585 -1.72116940 -1.74236542 -0.5169085 -0.5585368 -0.9516725 -1.205271 0.8613053 0.8425803 0.6477345 0.6692036
No -1.047396 -1.025007 -0.03288072 0.03602869 -1.0043711 -1.0653740 -0.9516725 -1.226141 0.7644880 0.7859356 0.5836250 0.5843434
No -1.062811 -1.055877 -0.81568941 -0.36312211 -1.2977547 -1.3621731 -0.9516725 -1.219184 0.8209648 0.8021198 0.7118441 0.7116337

Los datos han sido preparados y transformados, se continua con la selección de variables

Seleccion de variables

Se realiza una selección automática de variables de tipo stepwise para escoger las features que son mas relevantes para el posterior proceso de modelamiento y tuneado de algoritmos.

# Seleccion de variables -----
full <- glm(factor(varObjBin) ~., data = base, family = binomial(link = "logit"))
null <- glm(factor(varObjBin) ~ 1, data = base, family = binomial(link = "logit"))

# Seleccion de variables automatica
seleccion <- stepAIC(null, scope = list(upper = full), direction = "both", trace = F)

# Para ver las features escogidas de forma automatica
dput(names(seleccion$coefficients))
c("(Intercept)", "h10pix", "temp90", "Ymin", "Ymax", "temp", 
"humid", "humid90", "trees", "trees90")
# Version formula
formula(seleccion)
factor(varObjBin) ~ h10pix + temp90 + Ymin + Ymax + temp + humid + 
    humid90 + trees + trees90
summary(seleccion)

Call:
glm(formula = factor(varObjBin) ~ h10pix + temp90 + Ymin + Ymax + 
    temp + humid + humid90 + trees + trees90, family = binomial(link = "logit"), 
    data = base)

Deviance Residuals: 
     Min        1Q    Median        3Q       Max  
-2.63747  -0.15358  -0.02882   0.44958   3.04419  

Coefficients:
            Estimate Std. Error z value             Pr(>|z|)    
(Intercept)  -1.8848     0.1866 -10.102 < 0.0000000000000002 ***
h10pix        3.1524     0.2542  12.403 < 0.0000000000000002 ***
temp90        3.1225     0.7286   4.286            0.0000182 ***
Ymin         -1.8854     1.1527  -1.636             0.101914    
Ymax          1.6698     1.1480   1.455             0.145781    
temp         -2.6276     0.6948  -3.782             0.000156 ***
humid         3.7731     0.9470   3.984            0.0000677 ***
humid90      -3.4328     0.9786  -3.508             0.000452 ***
trees        -0.6970     0.1991  -3.500             0.000465 ***
trees90       0.5499     0.2221   2.477             0.013266 *  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 2695.4  on 1985  degrees of freedom
Residual deviance: 1033.8  on 1976  degrees of freedom
AIC: 1053.8

Number of Fisher Scoring iterations: 7

Según el código anterior, las variables mas relevantes seleccionadas son: h10pix + temp90 + Ymin + Ymax + temp + humid + humid90 + trees + trees90, con las cuales se aplica la función steprepetidobinaria() con criterio “AIC” para remuestreo:

# Seleccion AIC
listaAIC <- steprepetidobinaria(
    data = data,
    vardep = vardep,
    listconti = listconti,
    sinicio = 12345,
    sfinal = 12355,
    porcen = 0.8,
    criterio = "AIC")

tabla1 <- listaAIC[[1]]
knitr::kable(tabla1, "pipe")
modelo Freq contador
1 h10pix+temp90+trees+trees90+Ymin+Ymax+temp+humid+humid90 4 9
5 h10pix+temp90+trees+trees90+Ymin+temp+humid+humid90 3 8
8 h10pix+temp90+trees+trees90+temp+humid+humid90 2 7
10 h10pix+temp90+trees+trees90+Ymin+Ymax 1 6
11 h10pix+trees+Ymin+humid 1 4
dput(listaAIC[[2]][[1]])
c("h10pix", "temp90", "trees", "trees90", "Ymin", "Ymax", "temp", 
"humid", "humid90")
dput(listaAIC[[2]][[2]])
c("h10pix", "temp90", "Ymin", "temp", "humid", "humid90", "trees", 
"trees90")

Se repite el mismo proceso pero esta vez con criterio “BIC”:

# Seleccion con Bic
listaBIC <- steprepetidobinaria(
    data = data,
    vardep = vardep,
    listconti = listconti,
    sinicio = 12345,
    sfinal=12355, porcen = 0.8, criterio = "BIC")

tabla2 <- listaBIC[[1]]
knitr::kable(tabla2, "pipe")
modelo Freq contador
1 h10pix+temp90 8 2
9 h10pix+humid 1 2
10 h10pix+humid+trees 1 3
11 h10pix+temp90+Ymin 1 3
dput(listaBIC[[2]][[1]])
c("h10pix", "temp90")
dput(listaBIC[[2]][[2]])
c("h10pix", "humid")

Observando los resultados anteriores y obteniendo los set de variables se comparan con validación cruzada repetida:

# Seleccion de variables AIC
medias1 <- cruzadalogistica(
    data = data,
    vardep = "varObjBin", 
    listconti = c("h10pix", "temp90", "trees", "trees90", 
    "Ymin", "Ymax", "temp", "humid", "humid90"),
    listclass = c(""), 
    grupos = 4, sinicio = 1234, repe = 5)

medias1$modelo <- "Logística1"

# Seleccion de variables BIC
medias2 <- cruzadalogistica(
    data = data,
    vardep="varObjBin", 
    listconti = c("h10pix", "temp90"),
    listclass = c(""), 
    grupos = 4, 
    sinicio = 1234, 
    repe = 5)

medias2$modelo <- "Logística2"

La función genera_graficos() toma las medias de los modelos generados y entrega el respectivo gráfico de tipo boxplot para facilitar la visualización (se agrega como anexo):

genera_graficos(medias1, medias2)

Observaciones:

Algoritmos de ML para clasificación binaria

Redes neuronales

Tunning de redes neuronales

Una vez seleccionada las variables utilizando StepwiseAIC para logística, se comienza con el tunning de redes neuronales para un problema de clasificación binaria con la finalidad de encontrar los parámetros que se ajusten mejor para el dataset con el que se esta trabajando, realizar un nuevo modelamiento y comparar en este caso con regresión logística y mas adelante, con los próximos algoritmos a modelar (Bagging, Random Forest etc).

Tunning de redes neuronales con base accuracy

Utilizando las variables seleccionadas bajo StepwiseAIC se prueba un primer modelo:

set.seed(12345)
control_redes <- trainControl(method = "repeatedcv", number=4, repeats=5,
    savePredictions = "all")

avnnetgrid <- expand.grid(
    size = c(5, 10, 15, 20),
    decay = c(0.01, 0.1, 0.001), bag = FALSE)

# Con las variables seleccionadas stepwiseAIC
redavnnet1 <- train(
    varObjBin ~ h10pix + temp90 + trees + trees90 + Ymin + Ymax + temp + 
    humid + humid90,
    data = dengue,
    method = "avNNet", 
    linout = FALSE, 
    maxit = 100,
    trControl = control_redes, 
    tuneGrid = avnnetgrid,
    repeats = 5,
    trace = F)

knitr::kable(redavnnet1$results, "pipe")

Con el código anterior se ha generado la siguiente tabla la cual puede ser de utilidad para determinar cuales parámetros podrían ser útiles al momento de elegir para entrenar el modelo:

size decay bag Accuracy Kappa AccuracySD KappaSD
5 0.001 FALSE 0.9052328 0.8075327 0.0139857 0.0278779
5 0.010 FALSE 0.9229599 0.8423037 0.0131987 0.0269257
5 0.100 FALSE 0.9179245 0.8325153 0.0104558 0.0211637
10 0.001 FALSE 0.9207424 0.8384348 0.0132498 0.0265557
10 0.010 FALSE 0.9260815 0.8484666 0.0098502 0.0203012
10 0.100 FALSE 0.9197360 0.8361295 0.0099104 0.0201965
15 0.001 FALSE 0.9211444 0.8389624 0.0114062 0.0230180
15 0.010 FALSE 0.9266882 0.8494893 0.0102616 0.0213015
15 0.100 FALSE 0.9225561 0.8418606 0.0098902 0.0201859
20 0.001 FALSE 0.9244712 0.8454846 0.0101698 0.0207285
20 0.010 FALSE 0.9274934 0.8511611 0.0100521 0.0207691
20 0.100 FALSE 0.9216511 0.8400228 0.0101441 0.0207448

Observando la tabla y con base el Accuracy, los parámetros a considerar son un size de 10 con un decay de 0.010 o

Se prueba el modelo otra vez pero con las variables seleccionadas con criterio BIC (son solo dos las predictoras mas la variable objetivo)

# Con las variables seleccionadas BIC
redavnnet2 <- train(
    varObjBin ~ h10pix + temp90 ,
    data = data,
    method = "avNNet", linout = FALSE, maxit = 100,
    trControl = control_redes, tuneGrid = avnnetgrid,
    repeats = 5,
    trace = F)

knitr::kable(redavnnet2$results, "pipe")

La siguiente tabla muestra los resultados con las variables obtenidas con BIC

size decay bag Accuracy Kappa AccuracySD KappaSD
5 0.001 FALSE 0.9052328 0.8075327 0.0139857 0.0278779
5 0.010 FALSE 0.9229599 0.8423037 0.0131987 0.0269257
5 0.100 FALSE 0.9179245 0.8325153 0.0104558 0.0211637
10 0.001 FALSE 0.9207424 0.8384348 0.0132498 0.0265557
10 0.010 FALSE 0.9260815 0.8484666 0.0098502 0.0203012
10 0.100 FALSE 0.9197360 0.8361295 0.0099104 0.0201965
15 0.001 FALSE 0.9211444 0.8389624 0.0114062 0.0230180
15 0.010 FALSE 0.9266882 0.8494893 0.0102616 0.0213015
15 0.100 FALSE 0.9225561 0.8418606 0.0098902 0.0201859
20 0.001 FALSE 0.9244712 0.8454846 0.0101698 0.0207285
20 0.010 FALSE 0.9274934 0.8511611 0.0100521 0.0207691
20 0.100 FALSE 0.9216511 0.8400228 0.0101441 0.0207448

Agregar observaciones

Tunning de redes neuronales con MAXIT

Otro criterio para realizar tunning de redes neuronales es utilizando un bucle for() que recorre un vector llamado listaiter que contiene las iteraciones a evaluar, buscando nuevamente los parámetros que podrían ser mas adecuados para este set de datos

listconti = c("h10pix", "temp90", "trees", "trees90", "Ymin", "Ymax", 
    "temp", "humid", "humid90")

data2 <- dengue[,c(listconti, vardep)]

control_maxit <- trainControl(method = "cv", number = 4, savePredictions = "all") 

set.seed(12345)
nnetgrid <-  expand.grid(size = c(5, 10), decay = c(0.01, 0.1, 0.001), bag = F)

completo <-data.frame()
listaiter<-c(50, 100, 200, 500, 1000, 2000, 3000)

for( iter in listaiter) 
{
    rednnet<- train(
        varObjBin~.,
        data = data2,
        method = "avNNet",
        linout = FALSE,
        maxit = iter,
        trControl = control_redes, 
        repeats = 5, 
        tuneGrid = nnetgrid,
        trace = F)
    # Añado la columna del parametro de iteraciones
    rednnet$results$itera <- iter
    # Voy incorporando los resultados a completo
    completo <- rbind(completo, rednnet$results)
}

completo <- completo[order(completo$Accuracy),]

ggplot(completo, aes(x = factor(itera), y = Accuracy, 
    color = factor(decay), pch = factor(size))) +
    theme_minimal() + 
    geom_point(position = position_dodge(width = 0.5), size = 3)

Observando el gráfico con los resultados del tunning ajustando maxit, los mejores resultados en accuracy los entrega con un size de 10. En cuanto al decay muestra mejor comportamiento con decay 0.001 y 0.1. Con las iteraciones los mejores resultados los muestra entre 500 y 3000 aunque por simpleza con 500 iteraciones esta bien.

Validación cruzada repetida para redes neuronales

A continuación y una vez obtenidos los parámetros mas adecuados utilizando tunning con criterio accuracy y controlando maxit = se generan las correspondientes validaciones cruzadas repetidas utilizando la función cruzadaavnnetbin()

# Validacion cruzada repetida con criterio accuracy
medias3 <-cruzadaavnnetbin(
    data = dengue,
    vardep = "varObjBin",
    listconti = c("h10pix", "temp90", "trees", "trees90", "Ymin", "Ymax", 
    "temp", "humid", "humid90"),
    listclass = c(""),
    grupos = 4,
    sinicio = 1234,
    repe = 5,
    size = c(10), 
    decay = c(0.010), 
    repeticiones = 5, 
    itera = 200)

medias3$modelo <- "Avnnet1"

medias4 <- cruzadaavnnetbin(
    data = dengue,
    vardep = "varObjBin", 
    listconti = c("h10pix", "temp90"),
    listclass = c(""),
    grupos = 4,
    sinicio = 1234, 
    repe = 5,
    size = c(15),
    decay = c(0.10),
    repeticiones = 5,
    itera = 200)

medias4$modelo <- "Avnnet2"

# Validacion cruzada repetida tunning maxit
medias5 <- cruzadaavnnetbin(
    data = data2,
    vardep = "varObjBin",
    listconti = c("h10pix", "temp90", "trees", "trees90", "Ymin", "Ymax", 
    "temp", "humid", "humid90"),
    listclass = c(""), 
    grupos = 4, 
    sinicio = 1234,
    repe = 5, 
    repeticiones = 5,
    itera = 1000,
    size = c(10),
    decay = c(0.01))

medias5$modelo <- "Red con maxit"

Comparación de medias y boxplot

Utilizando la función genera_gráficos() se generan los gráficos de tipo boxplot para comparar los modelos creados.

genera_graficos(medias1, medias2, medias3, medias4, medias5)

Observaciones:

  • Si se compara con los modelos generados bajo regresión logística, hay un aumento considerable del criterio AUC con redes neuronales, específicamente la que se obtuvo con las variables seleccionadas con criterio AIC mas el modelo que utilizo los parámetros que entrego el tunning con maxit.

  • De todos modos esta es una comparativa inicial y aun se deben entrenar nuevos algoritmos y modelados para tomar una decisión.

Bagging

Tunning de bagging

En este apartado se probaran algunas técnicas para realizar tunning de bagging, partiendo de un modelo “básico” utilizando las variables previamente seleccionadas con stepwise y al final de este apartado realizando pruebas de tunning sobre el tamaño muestral del dataset dengue.

Primer modelo de bagging con variables previamente seleccionadas

rfgrid <- expand.grid(mtry = c(6))
set.seed(12345)
control_bagging <- trainControl(
    method = "cv", 
    number = 4, 
    savePredictions = "all", 
    classProbs = TRUE)

bagging <- train(
    data = dengue, 
    factor(varObjBin) ~ h10pix + temp90 + trees + trees90 + Ymin + Ymax + temp +
    humid + humid90,
    method = "rf", 
    trControl = control_bagging, 
    tuneGrid = rfgrid, 
    linout = FALSE, 
    ntree=5000,
    sampsize = 200, 
    nodesize = 10, 
    replace = TRUE)

bagging$modelType
mtry Accuracy Kappa AccuracySD KappaSD
6 0.9007989 0.7987109 0.008227 0.0164896

A continuación se explora el Out of bag error para este modelo, buscando parámetros para optimizacion, tal como lo indican los apuntes, se debe utilizar la función de randonForest() para extraer los datos y presentar gráfico:

# Plotear
baggingbis <- randomForest(
    factor(varObjBin) ~ h10pix + temp90 + trees + trees90 + Ymin + Ymax + temp +
    humid + humid90, 
    data = dengue, 
    mtry = 6, 
    ntree = 5000, 
    sampsize = 300, 
    nodesize = 10, 
    replace = TRUE)

#Se plotea la columna 1 que es OOB
plot(baggingbis$err.rate[, 1], main = "Out of bag error", ylab = "error")



# 1.3 Validacion cruzada rep bagging ------------------------------------------

A partir de las 1000 muestras se estabiliza el error de los datos analizados

Validacion cruzada repetida para bagging

Se realizan algunas pruebas con validación cruzada repetida utilizando la función cruzadarfbin() y cuidando de incluir el parámetro mtry

# Con validacion cruzada y 4 grupos
medias6 <-cruzadarfbin(
    data = dengue,
    vardep = "varObjBin",
    listconti = c("h10pix", "temp90", "trees", "trees90", "Ymin", "Ymax", 
    "temp", "humid", "humid90"),
    listclass = c(""), 
    grupos = 4, 
    sinicio = 1234, 
    repe = 5,
    mtry = 6, 
    ntree =1000, 
    replace = TRUE)

medias6$modelo <- "bagging"

A continuación y a modo de complemento se prueba a manipular el tamaño muestral para observar el efecto sobre bagging utilizando el parámetro sampsize. Para el dataset dengue son 2000 observaciones con 4 grupos de validación cruzada cada grupo contara con 1500 observaciones para probar:

Nota: El bagging medias7 o bagging_base es solo como adicional pues se hace con 10 grupos de validación cruzada, es solo para observar el efecto y comparar. medias8, medias9 y medias10 si se consideran pues utilizan 4 grupos de validación y con diversos tamaños muestrales:

#  Anexo Manipulando tamanio muestral con 10 grupos no incluir en la seleccion
#  de modelos pues el grupos = es diferente
medias7 <-cruzadarfbin(
    data = dengue,
    vardep = "varObjBin",
    listconti = c("h10pix", "temp90", "trees", "trees90", "Ymin", "Ymax", 
    "temp", "humid", "humid90"),
    listclass = c(""),
    grupos = 10,
    sinicio = 1234,
    repe=20, 
    nodesize = 10,
    mtry = 6, 
    ntree = 3000, 
    replace = TRUE)

medias7$modelo <- "bagging_base"

# Bagging manipulando sampsize 1000
medias8 <-cruzadarfbin(
    data = dengue,
    vardep = "varObjBin",
    listconti = c("h10pix", "temp90", "trees", "trees90", "Ymin", "Ymax", 
    "temp", "humid", "humid90"),
    listclass = c(""),
    grupos = 10,
    sinicio = 1234,
    repe=20, 
    nodesize = 10,
    mtry = 6, 
    ntree = 3000, 
    replace = TRUE, 
    sampsize = 1000 )

medias8$modelo <- "bagging1000"

# Bagging manipulando sampsize 1250
medias9 <-cruzadarfbin(
    data = dengue,
    vardep = "varObjBin",
    listconti = c("h10pix", "temp90", "trees", "trees90", "Ymin", "Ymax", 
    "temp", "humid", "humid90"),
    listclass = c(""),
    grupos = 10,
    sinicio = 1234,
    repe=20, 
    nodesize = 10,
    mtry = 6, 
    ntree = 3000, 
    replace = TRUE, 
    sampsize = 1250 )

medias9$modelo <- "bagging1250"

# Bagging manipulando sampsize 1500
medias10 <-cruzadarfbin(
    data = dengue,
    vardep = "varObjBin",
    listconti = c("h10pix", "temp90", "trees", "trees90", "Ymin", "Ymax", 
    "temp", "humid", "humid90"),
    listclass = c(""),
    grupos = 10,
    sinicio = 1234,
    repe=20, 
    nodesize = 10,
    mtry = 6, 
    ntree = 3000, 
    replace = TRUE, 
    sampsize = 1500 )

medias10$modelo <- "bagging1500"

Como parte de este complemento se evalúan las medias con la función genera_gráficos() para comparar cual de todos estos tamaños muestrales entrega mejores métricas:

# genera graficos para comparar tamanio muestral de bagging
genera_graficos(medias6, medias7, medias8, medias9, medias10, medias10)

Observaciones

  • De todo el grupo y descartando bagging_base que no se considera, los mejores resultados están entre bagging 1250 y bagging1500. Para comparar modelos en el apartado final incluyendo los modelos de random forest se deja el modelo medias9.

Random Forest Classifier

Tunning de random forest

Primer modelo utilizando todas las variables y buscando las variables predictoras mas relevantes.

Nota: Solo como complemento a la practica, para observar como se comporta el algoritmo al seleccionar variables. Para comparar medias, validación cruzada y probar con diferentes números de arboles se utilizan las variables seleccionadas con stepwise

# 2. Tuneo random forest -----------------------------------------------------
set.seed(12345)
rfgrid_rf <- expand.grid(mtry = c(3, 4, 5, 6, 7, 8, 9, 10, 11))

control_rf <- trainControl(
    method = "cv",
    number = 4,
    savePredictions = "all",
    classProbs = TRUE)

random_forest <- train(
    factor(varObjBin) ~.,
    data = dengue,
    method = "rf",
    trControl = control_rf,
    tuneGrid = rfgrid_rf,
    linout = FALSE,
    ntree =300,
    nodesize = 10,
    replace = TRUE,
    importance = TRUE)

# Mejor valor de mtry
random_forest$bestTune$mtry

Se obtienen los resultados de este primer modelo de randomforest

Para el modelo anterior el mejor mtry es 9 y en el gráfico se observa que las variables predictoras mas relevantes serian: Xmax, h10pix, Ymax, trees90, h10pix90, humid, Ymin lo cual difieren a los resultados entregados por la selección stepwise pero para explorar las alternativas que ofrece el algoritmo es útil.

A continuación se prueba generando un bucle for() para buscar los mejores indicadores evaluando según el mtry y numero de trees, desde 300 a 2500 para observar su comportamiento y buscar los mejores parámetros para la validación cruzada:

# random forest solo usando las variables mas relevantes
rfgrid_rf2 <- expand.grid(mtry = c(3, 4, 5, 6))
modellist <- list()
#train with different ntree parameters
for (ntree in c(300,500,1000,1500,2000,2500)){
    set.seed(12345)
    fit <- train(
        factor(varObjBin) ~ h10pix + temp90 + trees + trees90 + Ymin + Ymax + 
        temp + humid + humid90,
        data = dengue,
        method = "rf",
        trControl = control_rf,
        tuneGrid = rfgrid_rf2,
        lineout = FALSE,
        ntree = ntree,
        nodesize = 10,
        replace = TRUE,
        importance = TRUE)
    key <- toString(ntree)
    modellist[[key]] <- fit
}

results <- resamples(modellist)
dotplot(results)
dotplot(results)

fit$bestTune
  mtry
2    4

El gráfico dotplot() muestra el accuracy correspondiente a los trees evaluados en el ciclo for anterior, con esto los resultados en accuracy son similares pero por complejidad se prefiere dejar para pruebas en 500 trees. Para el parámetro bestTune el mtry es de 4 que también se dejara fijado en la validación cruzada.

Validacion cruzada para random forest

Se genera la validación cruzada para random forest utilizando la función cruzadarfbin y los parámetros obtenidos del estudio anterior, con mtry= 4 y ntree = 500

# Probar a dejar el parametro sampsize sin valor para que lo haga x defecto
medias11 <- cruzadarfbin(data = dengue, vardep = "varObjBin",
    listconti = c("h10pix", "temp90", "trees", "trees90", "Ymin", "Ymax",
    "temp", "humid", "humid90"),
    listclass = c(""),
    grupos = 4,
    sinicio = 1234,
    repe = 10,
    nodesize = 10,
    mtry = 4, ntree = 500, replace = TRUE)

medias11$modelo <- "randomforest"

Comparacion de medias y boxplot

Utilizando la función genera_gráficos() se generan los gráficos para realizar comparación de medias y evaluar los modelos

# 3. Compara medias ----------------------------------------------------------
genera_graficos(medias1, medias2, medias3, medias4, medias5,medias6, medias9,medias11)

Observaciones:

  • De los algoritmos evaluados hasta el momento, randomforest (que muestra mas varianza comparado con los otros algoritmos) y bagging son los que muestran mejor rendimiento en AUC y tasa de fallos, dejando atrás las regresiones logísticas y las redes neuronales.

Stochastic Gradient Boosting

Tunning de gradient boosting

Se realiza un primer tunning del algoritmo gradient boosting con caret, probando los parametros basicos de shrinkage, n.minobsinnode y n.trees mas las variables seleccionadas con stepwise

set.seed(12345)
# El shrinkage va desde los valores 0.0001 y 0.2 
# Cuanto mas alto mas rapido pero puede resultar menos fino
# n.minobsnode mide la complejidad
# interaction.depth = c(2) arboles binarios

gb_grid <- expand_grid(
    shrinkage = c(0.1, 0.05, 0.03, 0.01, 0.001),
    n.minobsinnode = c(5,10,20),
    n.trees = c(100, 500, 1000, 2000, 3000, 5000),
    interaction.depth = c(2))

control_gb <- trainControl(
    method = "cv", 
    number = 4, 
    savePredictions = "all", 
    classProbs = TRUE)

gbm <- train(
    factor(varObjBin) ~ h10pix + temp90 + trees + trees90 + Ymin + Ymax + temp +
    humid + humid90, 
    data = dengue, 
    method = "gbm", 
    trControl = control_gb,
    tuneGrid = gb_grid, 
    distribution = "bernoulli", 
    bag.fraction = 1, 
    verbose = FALSE)

gbm
plot(gbm)
gbm$bestTune
   n.trees interaction.depth shrinkage n.minobsinnode
51    1000                 2       0.1              5

Observaciones:

  • Segun gbm$bestTune los mejores resultados se consiguen con 1000 iteraciones, shrinkage de 0.1 y n.minobsinnode de 5. Los graficos muestran que en iteraciones desde 500 en adelante no se observan diferencias considerables y tienen a compartir resultados especialmente desde el shrinkage 0.06. El maximo accuracy estaria entre 0.06 y 0.10 de shrinkage ademas de minobsinnode 5 tal como lo indica bestTune.

Estudio de early stopping

Se prueba con algunos parametros para buscar en que iteracion se estabiliza el algoritmo

# Tuneado early stopping, se cambia bag.fraction por 0.5 en vez de 1
gbmgrid <- expand.grid(
    shrinkage = c(0.05), 
    n.minobsinnode = c(5),
    n.trees = c(50, 100, 300, 500, 800, 1000, 1200, 2000, 5000), 
    interaction.depth = c(2))

gbm2 <- train(
    factor(varObjBin) ~h10pix + temp90 + trees + trees90 + Ymin + Ymax + temp +
        humid + humid90,
    data = dengue, 
    method = "gbm", 
    trControl = control_gb,
    tuneGrid = gbmgrid, 
    distribution = "bernoulli", 
    bag.fraction = 0.5, 
    verbose = FALSE)

gbm2$bestTune
plot(gbm2, xlab = "Boosting iterations")

# Importancia variables
# summary(gbm2)
tabla <- summary(gbm2)
par(cex=1.5, las =2)
barplot(tabla$rel.inf, names.arg = row.names(tabla), col = "#F5E1FD", 
        main = "Importancia de variables - gradient boosting")
  n.trees interaction.depth shrinkage n.minobsinnode
9    5000                 2       0.1              5

Importancia de variables

Observaciones:

  • Pendiente, preguntar si se usan todas las variables o las seleccionadas con stepwise

Validacion cruzada para gradient boosting

Se genera la validación cruzada para gradient boosting utilizando los parametros obtenidos en el proceso anterior de tunning:

# 2.Validacion cruzada repetida gradient -----------------------------------------
# Segun el grafico entre 2000 y 5000 no hay tanta diferencia, se deja en 3000
medias12 <- cruzadagbmbin(
    data = dengue, 
    vardep="varObjBin", 
    listconti=c("h10pix", "temp90", "trees", "trees90", "Ymin", "Ymax",
    "temp", "humid", "humid90"), 
    listclass = c(""), 
    grupos = 4,
    sinicio = 1234,
    repe = 5, 
    n.minobsinnode = 5,
    shrinkage = 0.05,
    n.trees = 2000,
    interaction.depth = 2)

medias12$modelo <- "gbm"

Comparacion de medias y boxplot

Utilizando la función genera_gráficos() se generan los gráficos para realizar comparación de medias y evaluar los modelos

# 2.1 Compara medias -------
genera_graficos(medias1, medias2, medias3, medias4, medias5,medias6, medias11,
                medias12)

Observaciones:

  • gbm() tiene un buen auc y esta en el medio del grafico de la tasa de fallos, de todas maneras no alcanza los niveles que muestra un random forest o un bagging para este dataset.

Xgboost

Tunning de Xgboost

## 3. Tuneo Xgboost -----------
set.seed(12345) 
xgbmgrid <- expand.grid( 
    min_child_weight = c(5, 10, 20), 
    eta = c(0.1, 0.05, 0.03, 0.01, 0.001), 
    nrounds = c(100, 500, 1000, 5000), 
    max_depth = c(2, 3, 4, 5, 6), 
    gamma = 0, 
    colsample_bytree = 1, 
    subsample = 1)

control_xgboost <-trainControl(
    method = "cv", 
    number = 4, 
    savePredictions = "all", 
    classProbs = TRUE)

xgbm <- train(
    factor(varObjBin) ~h10pix + temp90 + trees + trees90 + Ymin + Ymax + temp +
    humid + humid90,
    data = dengue, 
    method = "xgbTree",
    trControl = control_xgboost,
    tuneGrid = xgbmgrid,
    verbose = FALSE, 
    verbosity = 0)

xgbm
plot(xgbm)
xgbm$bestTune

nrounds max_depth eta gamma colsample_bytree min_child_weight subsample
290 500 6 0.1 0 1 5 1

Observaciones:

  • Pendiente, es necesario hacer feature selection?

Estudio de early stopping con xgboost

# Early stopping
# Con los parametros entregados anteriormente en best tune
# Con las variables ya preseleccionadas
xgbmgrid2 <- expand.grid(
    eta = c(0.1), 
    min_child_weight = c(20), 
    nrounds = c(50, 100, 150, 200, 250, 300, 500, 1000), 
    max_depth = 6, 
    gamma = 0, 
    colsample_bytree = 1, 
    subsample = 1)

set.seed(12345) 
control_xgboost1 <-trainControl(
    method = "cv", 
    number = 4, 
    savePredictions = "all", 
    classProbs = TRUE)

xgbm2 <- train(
    factor(varObjBin) ~h10pix + temp90 + trees + trees90 + Ymin + Ymax + temp +
        humid + humid90,
    data = dengue, 
    method = "xgbTree", 
    trControl = control_xgboost1, 
    tuneGrid = xgbmgrid2, 
    verbose = FALSE)

plot(xgbm2)

# Se cambia la semilla para observar variaciones
set.seed(45673) 
control_xgboost2 <-trainControl(
    method = "cv", 
    number = 4, 
    savePredictions = "all", 
    classProbs = TRUE)

xgbm3 <- train(
    factor(varObjBin) ~h10pix + temp90 + trees + trees90 + Ymin + Ymax + temp +
        humid + humid90,
    data = dengue, 
    method = "xgbTree", 
    trControl = control_xgboost2, 
    tuneGrid = xgbmgrid2, 
    verbose=FALSE)

xgbm3$bestTune
plot(xgbm3)

# Importancia de variables
varImp(xgbm3)
plot(varImp(xgbm3))

nrounds max_depth eta gamma colsample_bytree min_child_weight subsample
7 500 6 0.1 0 1 5 1
nrounds max_depth eta gamma colsample_bytree min_child_weight subsample
7 500 6 0.1 0 1 5 1
# Solo con las variables seleccionadas en stepwise
xgbmgrid3 <- expand.grid(
    eta = c(0.1, 0.05, 0.03, 0.01, 0.001), 
    min_child_weight = c(5, 10, 20), 
    nrounds = c(50,100,150,200,250,300), 
    max_depth = c(2,3,4,5,6), 
    gamma = 0, 
    colsample_bytree = 1,
    subsample = 1)

set.seed(12345) 
xgbm4 <- train(
    factor(varObjBin) ~ h10pix + h10pix90 + Xmax + humid90 + 
    Ymax + temp90 + Xmin,
    data = dengue, 
    method = "xgbTree",
    trControl = control_xgboost, 
    tuneGrid = xgbmgrid3, 
    verbose = FALSE)

xgbm4$bestTune
plot(xgbm4)

Comparacion de medias y boxplot

Utilizando la función genera_gráficos() se generan los gráficos para realizar comparación de medias y evaluar los modelos

genera_graficos(medias1, medias2, medias3, medias4, medias5,medias6, medias11,
                medias12, medias13)

Observaciones:

  • Pendiente, es necesario hacer feature selection?

Support Vector Machines

Tunning de support vector machines o SVM

SVM Lineal

Para el SVM lineal se necesita la constante de regularización C. A continuación se construyen un par de modelos buscando un valor C optimo utilizando las variables seleccionadas con stepwise:

Como en los modelos anteriores, se utiliza el archivo dengue.RDS que contiene los datos preprocesados

# Carga de datos
dengue <- readRDS(file = "./data/processed/dengue.rds")
# 1. SVM -------------
# 1.1 SVM lineal -----

SVMgrid <- expand.grid(C = c(0.01, 0.05, 0.1, 0.2, 0.5, 1, 2, 5, 10))
set.seed(12345)
control_svm <- trainControl(method = "cv", number = 4, savePredictions = "all")

svm_1 <-train(
    data = dengue,
    factor(varObjBin) ~ h10pix + temp90 + trees + trees90 + Ymin + 
    Ymax + temp + humid + humid90, 
    method = "svmLinear", 
    trControl = control_svm, 
    tuneGrid = SVMgrid, 
    verbose = FALSE) 

Se muestran los resultados y el gráfico correspondiente

#svm_1
knitr::kable(svm_1$results, "pipe")
C Accuracy Kappa AccuracySD KappaSD
0.01 0.8826781 0.7639678 0.0019727 0.0047681
0.05 0.8907305 0.7804037 0.0061257 0.0123735
0.10 0.8907295 0.7801867 0.0073525 0.0145912
0.20 0.8917365 0.7820591 0.0071642 0.0142584
0.50 0.8937506 0.7860458 0.0069711 0.0137572
1.00 0.8947566 0.7880580 0.0082336 0.0165615
2.00 0.8942526 0.7869114 0.0090433 0.0185537
5.00 0.8932476 0.7848042 0.0076308 0.0159467
10.00 0.8932466 0.7847956 0.0084889 0.0177490
svm_1$bestTune
  C
6 1
plot(svm_1$results$C, svm_1$results$Accuracy, xlab = "C value", ylab = "Accuracy", 
     main = "Resultados svm lineal 1")

Observando bestTune y el gráfico, para este primer modelado de svm el valor C mas adecuado estaría entre 0.20 y 2 (donde se observa mejor accuracy en tabla), caret en bestTune indica un valor de C de 1 que es también lo que se observa en el gráfico. Con este nuevo valor de C se construye otro modelo donde se busca si es que hay un mejor valor de C en un grid de 0.1 hasta 1

SVMgrid_1 <- expand.grid(C = c(0.1, 0.2, 0.3, 0.4, 0.6, 0,7, 0.8, 0.9, 1))

svm_2 <- train(
    data = dengue, 
    factor(varObjBin) ~ h10pix + temp90 + trees + trees90 + Ymin + 
    Ymax + temp + humid + humid90, 
    method = "svmLinear", 
    trControl = control_svm, 
    tuneGrid = SVMgrid_1, 
    verbose = FALSE) 

Se despliegan los resultados

C Accuracy Kappa AccuracySD KappaSD
0.0 NaN NaN NA NA
0.1 0.8897386 0.7781448 0.0187763 0.0370724
0.2 0.8922568 0.7828897 0.0185648 0.0366145
0.3 0.8932628 0.7849067 0.0181678 0.0358879
0.4 0.8927598 0.7838505 0.0187761 0.0371822
0.6 0.8922568 0.7829515 0.0178185 0.0351353
0.8 0.8922568 0.7829515 0.0178185 0.0351353
0.9 0.8932638 0.7849717 0.0181591 0.0359239
1.0 0.8937689 0.7860518 0.0196068 0.0387306
7.0 0.8902447 0.7785518 0.0197410 0.0394319
  C
9 1

Para este dataset con SVM lineal el valor de C mas adecuado seria C=1.

SVM Polinomial

En este segundo apartado se utiliza un kernel polinomial buscando el valor de C, el grado del polinomio y la escala (C, degree y scale respectivamente), se entrena el modelo:

# 1.2 SVM Polinomial ----
SVMgrid_p <- expand.grid(
    C = c(0.01, 0.05, 0.1, 0.2, 0.5, 1, 2, 5, 10),
    degree = c(2, 3),
    scale = c(0.1, 0.5, 1, 2, 5))

svm_3 <- train(
    data = dengue, 
    factor(varObjBin) ~ h10pix + temp90 + trees + trees90 + 
    Ymin + Ymax + temp + humid + humid90, 
    method = "svmPoly", 
    trControl = control_svm, 
    tuneGrid = SVMgrid_p, 
    verbose = FALSE)

Se despliegan los resultados

degree scale C
49 3 2 0.5

En los resultados con bestTune se sugiere un degree de 3, scale 2 y un c value de 0.5. En el gráfico los mejores resultados en accuracy los muestra con degree 3, para c value muestra buenos resultados de 0.5 en adelante.

Como los mejores resultados los da el degree 3, se muestra un nuevo gráfico con ese grado y realizar observaciones:

En general scale tiende a generar una curva hasta el valor C=0.5 y después descender para volver a alcanzar un accuracy alto en C=10. Como un C muy grande puede tender al sobreajuste para este dataset es mejor un C 0,5 o 1 y escala 2

SVM RBF

En este tercer apartado se agrega al c value el parámetro sigma de varianza-escala buscando el valor gamma/sigma mas adecuado para este modelo. Gamma: a mayor gamma se puede presentar sobreajuste.

# 1.3 SVM RBF --------
#Se agrega el parametro sigma
SVMgrid_rbf <- expand.grid(
    C = c(0.01, 0.05, 0.1, 0.2, 0.5, 1, 2, 5, 10, 30),
    sigma = c(0.01, 0.05, 0.1, 0.2, 0.5, 1, 2, 5, 10, 30)
)

svm_4 <- train(
    data = dengue, factor(varObjBin) ~ h10pix + temp90 + trees + trees90 + 
    Ymin + Ymax + temp + humid + humid90, 
    method = "svmRadial", 
    trControl = control_svm, 
    tuneGrid = SVMgrid_rbf, 
    verbose = FALSE)

Se despliegan los resultados

sigma C
95 0.5 30

  • A mayor sigma (10, 30) no se observa que los resultados mejoren en términos de accuracy.

  • bestTune recomienda un sigma de 0.5 y un C de 30, pero desde c 1 hasta c 30 los resultados se observan relativamente parejos por lo que se mantiene el dejar el valor c en 1 y un sigma de 0.5

Validacion cruzada para support vector machines

Con los valores obtenidos de svm lineal, polinomial y RBF se realiza la validación cruzada para su posterior ploteo y comparación de medias para este modelo versus los generados con los algoritmos anteriores

# 2. Validacion cruzada rep SVM -------------------------------------------
#cv para lineal
medias14 <- cruzadaSVMbin(
    data = dengue, 
    vardep = "varObjBin",
    listconti = c("h10pix", "temp90", "trees", "trees90", 
    "Ymin", "Ymax", "temp", "humid", "humid90"),
    listclass = c(""),
    grupos = 4,
    sinicio = 1234,
    repe = 5, 
    C = 1
)

medias14$modelo <- "svmLineal"

# cv para poli
medias15 <- cruzadaSVMbinPoly(
    data = dengue, 
    vardep = "varObjBin",
    listconti = c("h10pix", "temp90", "trees", "trees90", 
    "Ymin", "Ymax", "temp", "humid", "humid90"),
    listclass = c(""),
    grupos = 4,
    sinicio = 1234,
    repe = 5, 
    C = 0.5,
    degree = 3,
    scale = 2
)

medias15$modelo <- "svmPoly"

# cv RBF
medias16 <- cruzadaSVMbinRBF(
    data = dengue, 
    vardep = "varObjBin",
    listconti = c("h10pix", "temp90", "trees", "trees90", 
    "Ymin", "Ymax", "temp", "humid", "humid90"),
    listclass = c(""),
    grupos = 4,
    sinicio = 1234,
    repe = 5, 
    C = 1,
    sigma = 0.5
)

medias16$modelo <- "svmRBF"

Comparacion de medias y boxplot

Utilizando la función genera_graficos() se despliega el gráfico de comparación de medias

# 3. Compara medias ----
genera_graficos(medias1,medias2,medias3, medias4, medias5,medias6, medias11, medias12, 
    medias13, medias14,medias15, medias16)

Observaciones

  • De los tres modelos el que presenta mayor varianza es el modelo polinomial, los resultados en accuracy sobre todo en svmRBF son buenos pero sigue siendo hasta ahora un randomForest o un xgbm (solo observando el gráfico) los que mejores resultados entregan.

Ensamblado de algoritmos

Nota: El código completo de la practica de ensamblados se encuentra en el siguiente repositorio de GitHub, para no extender innecesariamente el reporte: https://github.com/TiaIvonne/MasterUCM-MachineLearningR/blob/master/5_Ensamblado.R

Preparacion del ensamblado con caret ensemble

Antes de comenzar el modelado, se deben tener decididos los parámetros para los algoritmos que serán ensamblados (proceso obtenido anteriormente en el apartado de tunning) para los siguientes algoritmos:

Redes neuronales Gradient boosting machine Random Forest Classifier Support Vector Machines (lineal, polinomial y radial)

Con los resultados obtenidos se construyen los grids y se realiza el modelamiento

# Semilla
set.seed(12345)
repeticiones <- 10

# Evaluar los modelos
stackControl <- trainControl(
    method = "repeatedcv", 
    number = 4, 
    repeats = repeticiones,
    savePredictions = TRUE,
    classProbs = TRUE)

# grid gbm
gbmGrid <- expand.grid(
    n.trees = c(2000), 
    interaction.depth = c(2), 
    shrinkage =c(0.1), 
    n.minobsinnode = c(20))

# grid random forest 
rfGrid <- expand.grid(mtry=c(3))

# grid svm lineal
svmlinGrid <- expand.grid(C=c(0.08))

# grid svm ply
svmPolyGrid <- expand.grid(C=c(0.03),degree=c(3),scale=c(2))

# grid svm radial
svmRadialGrid <- expand.grid(sigma = c(0.5), C = c(30))

# Modelado
set.seed(12345)
models <- caretList(varObjBin~., 
    data = dengue, 
    trControl = stackControl,
    tuneList = list( 
        parrf = caretModelSpec(method = "rf", maxnodes = 30, n.trees = 200, nodesize = 10, sampsize = 150, tuneGrid = rfGrid), 
        glm = caretModelSpec(method = "glm"), 
        gbm = caretModelSpec(method = "gbm", tuneGrid = gbmGrid), 
        svmlinear = caretModelSpec(method = "svmLinear", tuneGrid = svmlinGrid), 
        svmPoly = caretModelSpec(method = "svmPoly", tuneGrid = svmPolyGrid), 
        svmradial = caretModelSpec(method = "svmRadial", tuneGrid = svmRadialGrid)))

Se despliegan los resultados obtenidos para el ensamblado

results <- resamples(models)
summary(results)
dotplot(results)

modelCor(results)
splom(results)
results[[2]]

ensemble <- caretEnsemble(models)
# Aquí se recomiendan los pesos para el ensamblado # de todos los modelos y se 
# ve la tasa de aciertos de cada modelo y ensamblado

summary(ensemble)

# 2. Validacion cruzada y kit ensamblado ----
# 2.1 Leer funciones ------
# 2.2 Prepapar archivo, variables , semilla y repeticiones ----
dput(names(dengue)) 

set.seed(12345) 

archivo <- dengue 

vardep <- "varObjBin" 
listconti <- c("h10pix", "temp90", "trees", "trees90", "Ymin", "Ymax", 
    "temp", "humid", "humid90") 

listclass <- c("") 
grupos <- 4 
sinicio <- 1234 
repe <- 15


# 2.3 Obtener datos de cv repetida para cada algoritmo y 
# procesar resultado ----
medias_1<-cruzadalogistica(data = archivo,
    vardep = vardep,
    listconti=listconti,
    listclass = listclass, grupos = grupos,
    sinicio = sinicio, repe = repe)

medias1bis <- as.data.frame(medias_1[1])
medias1bis$modelo<-"logistica"
predi1<-as.data.frame(medias_1[2])
predi1$logi<-predi1$Yes

medias_2<-cruzadaavnnetbin(data=archivo,
    vardep=vardep,listconti=listconti,
    listclass=listclass,grupos=grupos,sinicio=sinicio,repe=repe,
    size=c(10),decay=c(0.01),repeticiones=5,itera=200)

medias2bis<-as.data.frame(medias_2[1])
medias2bis$modelo<-"avnnet"
predi2<-as.data.frame(medias_2[2])
predi2$avnnet<-predi2$Yes

medias_3<-cruzadarfbin(data=archivo,
    vardep=vardep,listconti=listconti,
    listclass=listclass,grupos=grupos,sinicio=sinicio,repe=repe,
    mtry=3,ntree=500,nodesize=10,replace=TRUE)

medias3bis<-as.data.frame(medias_3[1])
medias3bis$modelo<-"randomforest"
predi3<-as.data.frame(medias_3[2])
predi3$rf<-predi3$Yes

medias_4 <-cruzadagbmbin(data=archivo,
    vardep=vardep,listconti=listconti,
    listclass=listclass,grupos=grupos,sinicio=sinicio,repe=repe,
    n.minobsinnode=5,shrinkage=0.1,n.trees=3000,interaction.depth=2)

medias4bis<-as.data.frame(medias_4[1])
medias4bis$modelo<-"gbm"
predi4<-as.data.frame(medias_4[2])
predi4$gbm<-predi4$Yes

medias_5 <-cruzadaxgbmbin(data = archivo,
    vardep = vardep,listconti = listconti,
    listclass = listclass, grupos = grupos, sinicio = sinicio, repe = repe,
    min_child_weight = 5, eta = 0.10, nrounds = 250, max_depth = 6,
    gamma = 0, colsample_bytree = 1, subsample = 1,
    alpha = 0, lambda = 0, lambda_bias = 0)

medias5bis <-as.data.frame(medias_5[1])
medias5bis$modelo <-"xgbm"
predi5 <-as.data.frame(medias_5[2])
predi5$xgbm<-predi5$Yes

medias_6<-cruzadaSVMbin(data=archivo,
    vardep=vardep,listconti=listconti,
    listclass=listclass,grupos=grupos,
    sinicio=sinicio,repe=repe,C=0.08)

medias6bis<-as.data.frame(medias_6[1])
medias6bis$modelo<-"svmLinear"
predi6<-as.data.frame(medias_6[2])
predi6$svmLinear<-predi6$Yes

medias_7 <- cruzadaSVMbinPoly(data = archivo,
     vardep = vardep, listconti = listconti,
     listclass = listclass, grupos = grupos, sinicio = sinicio, repe = repe,
     C = 0.03, degree = 3, scale = 2)

medias7bis <- as.data.frame(medias_7[1])
medias7bis$modelo <- "svmPoly"
predi7 <- as.data.frame(medias_7[2])
predi7$svmPoly <- predi7$Yes

medias_8 <- cruzadaSVMbinRBF(data = archivo,
    vardep = vardep, listconti = listconti,
    listclass = listclass, grupos = grupos,
    sinicio = sinicio, repe = repe,
    C = 30, sigma = 0.5)

medias8bis<-as.data.frame(medias_8[1])
medias8bis$modelo<-"svmRadial"
predi8<-as.data.frame(medias_8[2])
predi8$svmRadial<-predi8$Yes
summary(ensemble)
The following models were ensembled: parrf, glm, gbm, svmlinear, svmPoly, svmradial 
They were weighted: 
4.1057 -0.3249 0.4592 -2.7462 -1.4311 -1.5155 -3.104
The resulting Accuracy is: 0.9468
The fit for each individual model on the Accuracy is: 
    method  Accuracy  AccuracySD
     parrf 0.8973276 0.013143839
       glm 0.8899252 0.012870074
       gbm 0.9378647 0.009223015
 svmlinear 0.8885664 0.011650506
   svmPoly 0.9256787 0.011685866
 svmradial 0.9385687 0.009063826

Observaciones

Validacion cruzada repetida y boxplot

Siguiendo la guia de ensamblados entregada en el material, se realizan los pasos requeridos para realizar pruebas de ensamblado y validacion cruzada repetida

Paso 1, carga del archivo con la funcion “cruzadas ensamblado binaria fuente.R”

Paso 2, Preparacion del archivo

Se define semilla, variables y repeticiones

Paso 3, aplicacion de la funcion cruzadas ensamblado

Con los datos obtenidos del proceso de tunning se construyen las nuevas medias a evaluar bajo la funcion “cruzadas ensamblado binaria fuente.R”.

Solo se adjunta una muestra del código generado

medias_1<-cruzadalogistica(data = archivo,
    vardep = vardep,
    listconti=listconti,
    listclass = listclass, grupos = grupos,
    sinicio = sinicio, repe = repe)

medias1bis <- as.data.frame(medias_1[1])
medias1bis$modelo<-"logistica"
predi1<-as.data.frame(medias_1[2])
predi1$logi<-predi1$Yes

medias_2<-cruzadaavnnetbin(data=archivo,
    vardep=vardep,listconti=listconti,
    listclass=listclass,grupos=grupos,sinicio=sinicio,repe=repe,
    size=c(10),decay=c(0.01),repeticiones=5,itera=200)

medias2bis<-as.data.frame(medias_2[1])
medias2bis$modelo<-"avnnet"
predi2<-as.data.frame(medias_2[2])
predi2$avnnet<-predi2$Yes

medias_3<-cruzadarfbin(data=archivo,
    vardep=vardep,listconti=listconti,
    listclass=listclass,grupos=grupos,sinicio=sinicio,repe=repe,
    mtry=3,ntree=500,nodesize=10,replace=TRUE)

medias3bis<-as.data.frame(medias_3[1])
medias3bis$modelo<-"randomforest"
predi3<-as.data.frame(medias_3[2])
predi3$rf<-predi3$Yes

# Se ha omitido el resto del codigo

Con la función genera graficos y una vez calculadas las medias en el proceso anterior se obtiene el grafico de cajas respectivo:

# Genera graficos
genera_graficos(medias1bis, medias2bis, medias3bis, medias4bis, medias5bis, 
                medias6bis, medias7bis, medias8bis)

Paso 4, construccion del ensamblado

Con las predicciones obtenidas en las respectivas variables predi1, predi2, predi3 etc se crean los ensamblados. Solo se adjunta una muestra del codigo

unipredi<-cbind(predi1,predi2,predi3,predi4,predi5,predi6,predi7,predi8)
ncol(unipredi)

unipredi<- unipredi[, !duplicated(colnames(unipredi))]
ncol(unipredi)

unipredi$predi9<-(unipredi$logi+unipredi$avnnet)/2
unipredi$predi10<-(unipredi$logi+unipredi$rf)/2
unipredi$predi11<-(unipredi$logi+unipredi$gbm)/2
unipredi$predi12<-(unipredi$logi+unipredi$xgbm)/2
unipredi$predi13<-(unipredi$logi+unipredi$svmLinear)/2
unipredi$predi14<-(unipredi$logi+unipredi$svmPoly)/2

Paso 5, Procesado de los ensamblados

Solo se adjunta una parte del codigo, se construyen los promedios de tasa de fallos y AUC.

dput(names(unipredi))
listado<-c("logi", "avnnet", 
           "rf","gbm",  "xgbm", "svmLinear",  "svmPoly", 
           "svmRadial","predi9", "predi10", "predi11", "predi12", 
           "predi13", "predi14", "predi15", "predi16", "predi17", "predi18", 
           "predi19", "predi20", "predi21", "predi22", "predi23", "predi24", 
           "predi25", "predi26", "predi27", "predi28", "predi29", "predi30", 
           "predi31", "predi32", "predi33", "predi34", "predi35", "predi36", 
           "predi37", "predi38", "predi39", "predi40", "predi41", "predi42", 
           "predi43", "predi44", "predi45", "predi46", "predi47", "predi48", 
           "predi49", "predi50", "predi51", "predi52", "predi53", "predi54", 
           "predi55", "predi56", "predi57", "predi58", "predi59", "predi60", 
           "predi61", "predi62", "predi63", "predi64", "predi65", "predi66", 
           "predi67", "predi68", "predi69")

# Cambio a Yes, No, todas las predicciones
# Defino funcion tasafallos

tasafallos<-function(x,y) {
    confu<-confusionMatrix(x,y)
    tasa<-confu[[3]][1]
    return(tasa)
}

auc<-function(x,y) {
    curvaroc<-roc(response=x,predictor=y)
    auc<-curvaroc$auc
    return(auc)
}

# Se obtiene el numero de repeticiones CV y se calculan las medias por repe en
# el data frame medias0

repeticiones<-nlevels(factor(unipredi$Rep))
unipredi$Rep<-as.factor(unipredi$Rep)
unipredi$Rep<-as.numeric(unipredi$Rep)

medias0<-data.frame(c())
for (prediccion in listado)
{
    unipredi$proba<-unipredi[,prediccion]
    unipredi[,prediccion]<-ifelse(unipredi[,prediccion]>0.5,"Yes","No")
    for (repe in 1:repeticiones)
    {
        paso <- unipredi[(unipredi$Rep==repe),]
        pre<-factor(paso[,prediccion])
        archi<-paso[,c("proba","obs")]
        archi<-archi[order(archi$proba),]
        obs<-paso[,c("obs")]
        tasa=1-tasafallos(pre,obs)
        t<-as.data.frame(tasa)
        t$modelo<-prediccion
        auc<-suppressMessages(auc(archi$obs,archi$proba))
        t$auc<-auc
        medias0<-rbind(medias0,t)
    }
}

Paso 6, boxplot inicial

Boxplot con tasa de fallos para todos los ensamblados, en el punto 8 se muestran ordenados

Paso 7, tabla con resultados

Se genera tabla con resultados para la tasa de fallos y el area under curve

tablamedias<-medias0 %>%
    group_by(modelo) %>%
    summarise(tasa=mean(tasa))
    
tablamedias<-as.data.frame(tablamedias[order(tablamedias$tasa),])
knitr::kable(head(tablamedias, n = 10), "pipe")
modelo tasa
predi55 0.0641826
predi16 0.0650554
rf 0.0652904
predi57 0.0656260
predi26 0.0663310
predi60 0.0663646
predi63 0.0663981
predi69 0.0667338
predi52 0.0670695
predi56 0.0678751
# Para AUC
tablamedias2<-medias0 %>%
    group_by(modelo) %>%
    summarise(auc=mean(auc))     

tablamedias2<-tablamedias2[order(-tablamedias2$auc),]
knitr::kable(head(tablamedias2, n=10), "pipe")
modelo auc
predi55 0.9789378
predi57 0.9786039
predi52 0.9782842
predi60 0.9782593
predi23 0.9780586
predi69 0.9779774
predi56 0.9778908
predi26 0.9778800
predi16 0.9777369
predi54 0.9775005

En el punto siguiente se grafica para una mejor visualizacion y toma de decision

Paso 8, Boxplot ordenados para comparacion de medias

Con los resultados obtenidos de los ensamblados se generan los boxplot ordenados por tasa de fallos y auc respectivamente:

knitr::include_graphics("./figure/ensamblado_medias0_tasa.png")

En tasa de fallos los modelos con menor tasa son predi55, predi16, randomforest

knitr::include_graphics("./figure/ensamblado_medias0_auc.png")

En AUC respectivamente los modelos con mayor accuracy son predi55, predi57 y predi52.

Como se observa en el grafico son demasiados modelos, en el apartado siguiente solo se grafican los mejores modelos de ensamblado y se comparan con los alfgoritmos sin ensamblar

Paso 9, comparacion de mejores modelos

Los mejores ensamblados son los modelos predi57, predi55, predi52 y predi 60, estos se muestran a continuacion comparandolos con los modelos originales sin emsamblar

Paso 10, revision a los mejores ensamblados

unipredi$predi55<-(unipredi$rf+unipredi$xgbm+unipredi$svmRadial)/3
unipredi$predi57<-(unipredi$rf+unipredi$avnnet+unipredi$xgbm)/3
unipredi$predi60<-(unipredi$rf+unipredi$avnnet+unipredi$svmRadial)/3
unipredi$predi52<-(unipredi$rf+unipredi$gbm+unipredi$svmRadial)/3

En los ensamblados el algoritmo común en los cuatro mejores ensamblados es random forest, el mejor ensamblado es el numero 55 que mezcla random forest, xgbm y svm radial.

Observaciones finales

  • De los modelos originales los mejores resultados observando el grafico los obtiene random forest, xgbm y avnnet, xgbm presenta mas varianza que los dos anteriores.
  • El accuracy alcanzado por los ensamblados es mas alto si se compara con los modelos originales sin ensamblar. Si solo fuese tomando como criterio el accuracy, predi57 que es un ensamblado de randomforest con xgbm y svmlineal seria el modelo ganador. Sin embargo no hay diferencias dramáticas desde randomforest hacia adelante y considerando sesgo-varianza, randomforest estaria mostrando mejores resultados que los ensamblados.
  • La tasa de fallos vs el auc muestra algunas incoherencias, en tasa de fallos el numero menor lo obtienen predi55 y predi16 pero en auc los ensamblados predi55 y predi57 respectivamente estan a la cabecera de los resultados.

Evaluacion del modelo ganador y conclusiones

Algoritmo a usar para este modelo

Matriz de confusion para el modelo escogido

Sensitividad, especificidad y precision

Tabla de parametros para regresion logistica

Puntos de corte

Contraste de hipotesis entre algunos modelos

Visualpred

Tabla resumen

Como punto final a esta practica se adjunta tabla resumen con auc y tasa de fallos para cada modelo evaluado

lista_medias <- rbind(medias1,medias2,medias3, medias4, medias5,medias6, medias11, medias12, medias13, medias14,medias15, medias16)

knitr::kable(lista_medias, "pipe")
tasa auc modelo
0.1082578 0.9451795 Logística1
0.1072508 0.9466813 Logística1
0.1077543 0.9476129 Logística1
0.1112790 0.9467398 Logística1
0.1062437 0.9478490 Logística1
0.1117825 0.9440390 Logística2
0.1122860 0.9440421 Logística2
0.1122860 0.9449173 Logística2
0.1117825 0.9448233 Logística2
0.1127895 0.9448474 Logística2
0.0765358 0.9727161 Avnnet1
0.0730111 0.9716842 Avnnet1
0.0679758 0.9705218 Avnnet1
0.0730111 0.9714733 Avnnet1
0.0765358 0.9722795 Avnnet1
0.1122860 0.9442092 Avnnet2
0.1137966 0.9441967 Avnnet2
0.1122860 0.9443377 Avnnet2
0.1153072 0.9447805 Avnnet2
0.1158107 0.9446468 Avnnet2
0.0720040 0.9729448 Red con maxit
0.0720040 0.9722033 Red con maxit
0.0694864 0.9720247 Red con maxit
0.0735146 0.9722242 Red con maxit
0.0750252 0.9718148 Red con maxit
0.0704935 0.9745558 bagging
0.0674723 0.9754180 bagging
0.0745217 0.9760263 bagging
0.0740181 0.9764947 bagging
0.0689829 0.9757741 bagging
0.0634441 0.9776666 randomforest
0.0624371 0.9774295 randomforest
0.0664653 0.9782164 randomforest
0.0664653 0.9782368 randomforest
0.0659617 0.9781428 randomforest
0.0669688 0.9764060 randomforest
0.0634441 0.9761339 randomforest
0.0639476 0.9796274 randomforest
0.0669688 0.9759428 randomforest
0.0594159 0.9783731 randomforest
0.0901309 0.9707088 gbm
0.0845921 0.9673991 gbm
0.0855992 0.9687620 gbm
0.0896274 0.9691035 gbm
0.0795569 0.9705260 gbm
0.0503525 0.9883863 xgbm
0.0513595 0.9863748 xgbm
0.0528701 0.9893158 xgbm
0.0473313 0.9898651 xgbm
0.0468278 0.9899038 xgbm
0.1122860 0.9452891 svmLinear
0.1132931 0.9454886 svmLinear
0.1112790 0.9459857 svmLinear
0.1122860 0.9462740 svmLinear
0.1107754 0.9460902 svmLinear
0.0906344 0.9607452 svmPoly
0.0861027 0.9614794 svmPoly
0.0936556 0.9577070 svmPoly
0.0901309 0.9567107 svmPoly
0.0861027 0.9611755 svmPoly
0.0760322 0.9670231 svmRBF
0.0795569 0.9624653 svmRBF
0.0750252 0.9648027 svmRBF
0.0780463 0.9621196 svmRBF
0.0775428 0.9641823 svmRBF

Anexos

Otros metodos de seleccion de variables

Como un complemento a la selección de variables de tipo Stepwise que se ha estudiado en el apartado de selección de variables se ha decidido estudiar otras alternativas de selección de features utilizando diversas técnicas que se detallan a continuación:

Utilizando la libreria party y random forest

Un método de selección de variables es utilizando el algoritmo de random forest para encontrar un set de predictores, para eso se utiliza la librería party y se obtienen los siguientes resultados:

library(party)
s1 <- cforest(varObjBin ~ . , 
    data = dengue2, 
    control = cforest_unbiased(mtry = 2, ntree = 50))
varimp(s1)
      humid     humid90        temp      temp90      h10pix    h10pix90 
0.019011819 0.024480176 0.007135328 0.010393915 0.044306715 0.036816130 
      trees     trees90        Xmin        Xmax        Ymin        Ymax 
0.004198915 0.005095990 0.011027785 0.009075744 0.028311803 0.024437635 

Agregar un comentario

Utilizando MARS: Multivariate Adaptive Regression Splines

Otro método de selección de variables es utilizando MARS, contenida en la librería earth para R

library(earth)
s2 <- earth(varObjBin ~ ., data=dengue2) # build model
evaluacion <- evimp (s2)
plot(evaluacion)

Utilizando la libreria Boruta

Boruta es otra librería para realizar feature selection de forma automática, permitiendo un acercamiento inicial rápido al dataset (con sus ventajas y desventajas):

library(Boruta)
boruta_model <- Boruta(varObjBin ~ ., data= dengue2, doTrace=2)  
relevancia <- names(boruta_model$finalDecision[boruta_model$finalDecision %in% c("Confirmed")])  
print(relevancia) 
 [1] "humid"    "humid90"  "temp"     "temp90"   "h10pix"   "h10pix90"
 [7] "trees"    "trees90"  "Xmin"     "Xmax"     "Ymin"     "Ymax"    
plot(boruta_model, las = 2, cex.axis=0.7, xlab="")

Boruta permite revisar si el set de variables a estudiar entran en el modelo o podrían entrar de forma tentativa o definitivamente no han sido considerados como variables predictoras relevantes con lo cual también se pueden probar otros modelos y combinaciones (ej confirmados mas algún tentativo etc)

Para este dataset en particular la selección no marca variables como tentativas o descartadas por lo cual correspondería observar el gráfico generado y realizar algún corte de tipo manual en las variables.

Utilizando RFE con caret

Eliminación recursiva de características o RFE en sus siglas en ingles, es un estimador que asigna pesos a las características del dataset que se quiere estudiar

set.seed(7)
# Cargar datos
data(dengue2)
# Se define la funcion de control
controlrfe <- rfeControl(functions=rfFuncs, method="cv", number=4)
# Algoritmo RFE
res_rfe <- rfe(dengue2[,2:12], dengue2[,1], sizes=c(2:12), rfeControl=controlrfe)
#Resultados
print(res_rfe)

Recursive feature selection

Outer resampling method: Cross-Validated (4 fold) 

Resampling performance over subset size:

 Variables   RMSE Rsquared     MAE   RMSESD RsquaredSD    MAESD Selected
         2 0.2559   0.7144 0.13033 0.080328   0.188079 0.080015         
         3 0.2103   0.8186 0.09381 0.008193   0.013435 0.001988         
         4 0.2093   0.8200 0.09198 0.006153   0.011944 0.003804         
         5 0.1927   0.8479 0.08635 0.009504   0.014010 0.004464         
         6 0.1915   0.8499 0.08543 0.008331   0.011868 0.003713        *
         7 0.1934   0.8473 0.08939 0.006407   0.009145 0.004466         
         8 0.1941   0.8465 0.09210 0.004874   0.007074 0.003773         
         9 0.1949   0.8450 0.09183 0.004130   0.007247 0.003856         
        10 0.1985   0.8395 0.09571 0.004161   0.006483 0.003471         
        11 0.2012   0.8353 0.09888 0.002943   0.007213 0.002814         

The top 5 variables (out of 6):
   Xmax, h10pix, Xmin, h10pix90, Ymin
# Features escogidas
predictors(res_rfe)
[1] "Xmax"     "h10pix"   "Xmin"     "h10pix90" "Ymin"     "temp"    
# graficos
# plot(results, type=c("g", "o"))

Con las variables que ha seleccionado el modelo de RFE mas los modelos anteriores se prueban con la función cruzadalogistica() y se genera el boxplot para comparar:

# Generado con libreria party
anexo1 <- cruzadalogistica(
    data = data,
    vardep = "varObjBin", listconti = c("humid", "humid90","temp", "temp90",
    "h10pix", "h10pix90"),
    listclass = c(""), 
    grupos = 4, sinicio = 1234, repe = 5)

anexo1$modelo <- "log-party"


# Generado con MARS
anexo2 <- cruzadalogistica(
    data = data,
    vardep = "varObjBin", 
    listconti = c("h10pix", "Xmax", "temp90", "Ymax", "Ymin", "temp", "Xmin",
                  "trees"),
    listclass = c(""), 
    grupos = 4, sinicio = 1234, repe = 5)

anexo2$modelo <- "log-mars"

# Generado con Boruta
anexo3 <- cruzadalogistica(
    data = data,
    vardep = "varObjBin", 
    listconti = c("h10pix", "Xmax", "Xmin", "h10pix90", "Ymax"),
    listclass = c(""), 
    grupos = 4, sinicio = 1234, repe = 5)

anexo3$modelo <- "log-boruta"


# Generado con RFE
anexo4 <- cruzadalogistica(
    data = data,
    vardep = "varObjBin", 
    listconti = c("Xmax", "h10pix","Xmin","h10pix90","Ymin","temp" ),
    listclass = c(""), 
    grupos = 4, sinicio = 1234, repe = 5)

anexo4$modelo <- "log-RFE"

Conclusiones:

  • A modo general la selección de variables realizada con stepwiseAIC (Logistica1), sigue siendo la que presenta mejor AUC y menor tasa de fallos, ademas de menor varianza por lo que se seguirá trabajando con ese modelo.

Arboles de clasificacion

Seleccion de variables

Como complemento al trabajo realizado anteriormente con los algoritmos de ML y su posterior evaluación, se ha añadido un apartado de practica modelando arboles de clasificación utilizando R y el dataset dengue

Primer modelado, modelo con todas las variables del dataset para realizar estudio de importancia de variables (que se hizo en el apartado de selección con Stepwise pero se prefiere agregarlo a la practica)

Preguntar si en arboles trabajo con las variables que da el modelo o me ajusto a lo que dio la selección de variables con stepwise

# Modelo completo con gini para variable dependiente categorica
arbol1 <- rpart(factor(varObjBin) ~ ., 
    data = dengue_arbol, 
    minbucket = 30, 
    method = "class",
    parms = list(split = "gini"))

# Reglas de decision
rattle::asRules(arbol1)

Representación gráfica de las variables consideradas relevantes y árbol:

# Importancia de variables
par(cex = 0.7)
barplot(height = arbol1$variable.importance, col = "#e7e1fd", 
    main = "Importancia de variables Arbol1")

# Ploteo del arbol
rpart.plot(arbol1, extra = 105, tweak = 1.2, type = 1, nn = TRUE)

Observaciones:

  • Según el gráfico importancia de variables Árbol1, las varibles predictoras mas relevantes serian, h10pix, h10pix90, humid90, Ymin, humid e Ymax.

  • agregar observaciones del árbol

Se prueba modelando un nuevo árbol esta vez usando maxsurrogate = 0, solo como complemento puesto que el archivo dengue no contiene datos faltantes:

# Modelo con maxsurrogate = 0 para que trabaje los na's
arbol2 <- rpart(factor(varObjBin) ~ ., 
    data = dengue_arbol, 
    minbucket = 30, 
    method = "class", 
    maxsurrogate = 0,
    parms = list(split = "gini"))

Representación gráfica de las variables consideradas relevantes y árbol:

# Importancia de variables
barplot(height=arbol2$variable.importance, col = "#e7e1fd", 
        main = "Importancia de variables Arbol2")

# Ploteo del arbol
rpart.plot(arbol2, extra = 105, tweak = 1.2, type = 1, nn = TRUE)

Observaciones:

  • Según el gráfico importancia de variables Arbol2, las varibles predictoras mas relevantes serian, h10pix e Ymax.

  • agregar observaciones del árbol

Tunning de algoritmos usando rpart()

De forma inicial se realiza un primer tunning de arboles observando el comporta miento de la complejidad del árbol graficando su estructura en relación al parámetro minbucket (5, 30 y 60) utilizando la función rpart()

# Tuneado sobre complejidad del arbol con Rpart

arbol_min5 <- rpart(factor(varObjBin) ~ ., data = dengue_arbol, minbucket = 5, cp = 0)
arbol_min30 <- rpart(factor(varObjBin) ~ ., data = dengue_arbol, minbucket = 30, cp = 0)
arbol_min60 <- rpart(factor(varObjBin) ~ ., data = dengue_arbol, minbucket = 60, cp = 0)

Representación gráfica del tunning de minbucket

par(mfrow= c(2,2))
barplot(height = arbol_min5$variable.importance, col = "#e7e1fd", 
        main = "Importancia de variables con minbucket 5")
barplot(height = arbol_min30$variable.importance, col = "#e7e1fd", 
        main = "Importancia de variables con minbucket 30")
barplot(height = arbol_min60$variable.importance, col = "#e7e1fd", 
        main = "Importancia de variables con minbucket 60")
par(mfrow= c(2,2))

rpart.plot(arbol_min5, extra = 1)
rpart.plot(arbol_min30, extra = 1)
rpart.plot(arbol_min60, extra = 1)

Observaciones

Se puede apreciar que con minbucket 5 se genera un modelo complejo y con minbucket 60 un modelo básico, para encontrar un equilibrio y ver otras opciones de minbucket a continuación se profundiza en el tuneado utilizando caret de R

Tunning de algoritmos usando caret()

Como indican los apuntes, no se puede realizar tunning en el grid por lo que se construye un bucle for() que recorre por cada numero de minbucket indicado y va guardando los resultados para una posterior evaluación:

# En los ejemplos anteriores se dejaron los na's para ver el comportamiento, se eliminan 
# tunning de minubucket en bucle
# Tratamiento de datos faltantes
dengue_arbol1 <- na.omit(dengue_arbol, (!is.na(dengue_arbol)))
colSums(is.na(dengue_arbol1)) > 0 

# Tunning con minbucket en bucle
set.seed(12345)
control_arbol <- trainControl(method = "cv",
    number = 4,
    classProbs = TRUE, 
    savePredictions = "all") 

arbolgrid <-  expand.grid(cp = c(0))
tabla_resultados <-c() 

for (minbu in seq(from = 5, to = 60, by = 5)){
    arbolgrid <-  expand.grid(cp=c(0))
    arbolcaret <- train(factor(varObjBin) ~ h10pix + temp90 + Ymin + Ymax + 
        temp + humid + humid90 + trees + trees90,
    data = dengue_arbol1, 
    method = "rpart", 
    minbucket = minbu, 
    trControl = control_arbol, 
    tuneGrid = arbolgrid)
    
    accuracy <- arbolcaret$results$Accuracy
    sal <- arbolcaret$pred
    salconfu <- confusionMatrix(sal$pred, sal$obs)
    curvaroc <-roc(response = sal$obs, predictor = sal$Yes)
    auc <- curvaroc$auc
    # Guarda los resultados en formato tabla
    tabla_resultados <- rbind(tabla_resultados, c(minbu, accuracy, auc))
}

# Cambia los nombres de las columnas
colnames(tabla_resultados) <- c("minibucket","accuracy","AUC")

Con el código anterior se genera una tabla para una mejor visualización de los indicadores:

knitr::kable(tabla_resultados, "pipe")
minibucket accuracy AUC
5 0.9103816 0.9553222
10 0.8977636 0.9441889
15 0.8932587 0.9576898
20 0.9113887 0.9573321
25 0.9083675 0.9552511
30 0.8917355 0.9550814
35 0.8947678 0.9513164
40 0.9058403 0.9512422
45 0.9023131 0.9581003
50 0.9083574 0.9532965
55 0.9048383 0.9542511
60 0.9023232 0.9540407

Como se puede apreciar en la tabla anterior, en combinación entre mejor accuracy y AUC para este modelo es adecuado utilizar un minbucket de 30. Con ese numero se realiza la validación cruzada repetida utilizando la función cruzadaarbolbin()

arbol4 <- train(factor(varObjBin) ~ h10pix + temp90 + Ymin + Ymax + temp + 
    humid + humid90 + trees + trees90,
    data = dengue_arbol1, 
    method = "rpart", 
    minbucket = 20, 
    trControl = control_arbol, 
    tuneGrid = arbolgrid)

Se obtiene la matriz de confusión para este modelo

salida_a4 <- arbol4$pred
salida4_confusion <- confusionMatrix(salida_a4$pred, salida_a4$obs)
salida4_confusion
Confusion Matrix and Statistics

          Reference
Prediction   No  Yes
       No  1039   84
       Yes  123  740
                                               
               Accuracy : 0.8958               
                 95% CI : (0.8815, 0.9089)     
    No Information Rate : 0.5851               
    P-Value [Acc > NIR] : < 0.00000000000000022
                                               
                  Kappa : 0.7868               
                                               
 Mcnemar's Test P-Value : 0.008262             
                                               
            Sensitivity : 0.8941               
            Specificity : 0.8981               
         Pos Pred Value : 0.9252               
         Neg Pred Value : 0.8575               
             Prevalence : 0.5851               
         Detection Rate : 0.5232               
   Detection Prevalence : 0.5655               
      Balanced Accuracy : 0.8961               
                                               
       'Positive' Class : No                   
                                               
fourfoldplot(salida4_confusion$table)

Se obtiene la curva ROC

# Curva roc
curvaroc <-roc(response = salida_a4$obs, predictor = salida_a4$Yes)
auc <- curvaroc$auc
plot.roc(roc(response = salida_a4$obs, predictor = salida_a4$Yes), main = "Curva ROC arbol con minbucket = 20")

Comparacion con otros algoritmos de ML

Con los resultados obtenidos de entrenar el modelo se realiza la validación cruzada repetida para posteriormente comparar las medias con los demás modelos y evaluar sesgo-varianza

# 2. Validacion cruzada repetida -----
# Con los resultados de la tabla con iteraciones se configura la validacion cruzada
medias_arbol <-cruzadaarbolbin(
    data = dengue_arbol1,
    vardep = "varObjBin",
    listconti = c("h10pix", "temp90", "trees", "trees90", "Ymin", "Ymax", "temp", 
    "humid", "humid90"),
    listclass = c(""), grupos = 4, sinicio = 1234, repe = 5,
    cp = c(0), minbucket = 20)

medias_arbol$modelo <- "Arbol"

Con la función genera_gráficos() se despliega la comparación de medias para los distintos algoritmos entrenados

#genera_graficos(medias1, medias2, medias3, medias4, medias5,medias6, medias9,medias11, medias_arbol)

genera_graficos(medias1, medias2, medias_arbol)

Observaciones:

  • Para esta practica en particular, un árbol no es lo suficientemente competitivo frente a otros algoritmos modelados. Presenta mejores indicadores que la regresión logística o una red neuronal pero ante otros algoritmos como bagging o random forest se queda un tanto atrás.

Libreria h2o para python

Como practica adicional se ha decidido realizar un modelado con autoML utilizando la librería h2o para python, utilizando un jupyter notebook para este fin.

El notebook se puede revisar en detalle en: https://github.com/TiaIvonne/MasterUCM-MachineLearningR/blob/master/h2o%20python/autoML_dengue.ipynb

El primer paso ha sido generar un modelo de autoML que determine cual combinacion es la mas adecuada para el set de datos dengue.

En segundo lugar y con el modelo ganador escogido anteriormente, se hace un training del modelo para comparar resultados con los obtenidos en R

Codigo funcion genera_graficos()

La siguiente función toma como base lo visto en los apuntes con el fin de poder ser reutilizada y no estar redundando código en cada comparación de medias para cada evaluación de algoritmos:

genera_graficos <- function(...) {
    # Esta funcion toma las medias generadas con validacion cruzada, las ordena
    # y presenta dos graficos, segun tasa y auc
    
    # bind 
    union <-rbind(...)
    # ordena las medias segun tasa y auc
    ordenado <- union
    ordenado$modelo <- with(ordenado, reorder(modelo, tasa, median))
    ordenado2 <- union
    ordenado2$modelo <- with(ordenado2, reorder(modelo, auc, median))
    
    # colores para el grafico
    nb_cols <- 18
    mycolors <- colorRampPalette(brewer.pal(8, "Set2"))(nb_cols)
    
    # genera los graficos
    plot1 <- ggplot(ordenado, aes(x=modelo, y=tasa, fill = modelo)) + 
        ggtitle ("Tasa fallos")+
        geom_boxplot(alpha=0.3) +
        theme_minimal() + 
        theme(legend.position="none") +
        scale_fill_manual(values = mycolors) +
        xlab("Modelos") + ylab("Tasa") +
        theme(axis.text.x = element_text(angle = 45, hjust=1)) 
    
    plot2 <- ggplot(ordenado2, aes(x=modelo, y=auc, fill = modelo)) + 
        ggtitle ("Area Under Curve - AUC")+
        geom_boxplot(alpha=0.3) +
        theme_minimal() + 
        theme(legend.position="none") +
        scale_fill_manual(values = mycolors) +
        xlab("Modelos") + ylab("AUC") +
        theme(axis.text.x = element_text(angle = 45, hjust=1)) 
    list(plot1, plot2)
    # a. es donde se guarda la lista con plot1 y plot2 si se quiere
    # generar con grid.arrange
    # grid.arrange(grobs=a,ncol=2)
    
}

}